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Abstract 



o 

Cn In this paper, we provide a mathematical model for the electrolocation in weakly electric 

Jh fishes. We first investigate the forward complex conductivity problem and derive the approx- 

imate boundary conditions on the skin of the fish. Then we provide a dipole approximation 
for small targets away from the fish. Based on this approximation, we obtain a non-iterative 
location search algorithm using multi-frequency measurements. We present numerical experi- 
ments to illustrate the performance and the stability of the proposed multi-frequency location 
search algorithm. Finally, in the case of disk- and ellipse-shaped targets, we provide a method 

Q, to reconstruct separately the conductivity, the permittivity, and the size of the targets from 

multi-frequency measurements. 
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1 Introduction 



> 

QO In the turbid rivers of Africa and South America, some species of fish generate an electric current 

Cn which is not enough for defense purpose. In 1958, Lissmann and Machin discover that this electric 

^^ current is in fact used for spatial visualization "H^f . Indeed an object in the vicinity of the fish will 

be detected by measurement of the electric field's distortion on the skin. Behavioral experiments 
m have shown that the weakly electric fish is able to extract useful information about targets, such as 

^—^ the location fSO", the shape [35^, and the electric parameters (capacitance and conductivity) [i5] . 

Mathematically speaking, this is an inverse problem for the electric field created by the fish. 
Indeed, given the current distribution over the skin, the problem is to recover the conductivity 
distribution in the surrounding space. Due to the ill-posedness of this type of problems, it is very 
difficult to recover as much information as the fish is able to. Thus, modelling this "electric sense" 
$-H (called active electrolocation) is likely to give us insights in this regard. 

Electrolocation has been quantitatively investigated since Lissmann and Machin, who tried an 
analytical approach. More precisely, they computed the distortion created by a cylinder placed in 
the electric field of a dipole [32] , and noticed that it is equivalent to the field created by a dipole 
located inside the cylinder. In 1983, Bacher remarked that this formula cannot explain the phase 
difference observed when the electric permittivity of the target does not equal the permittivity of 
the water [13] . This phase shift seems to be an important input for the fish since it is measured by 
receptors (called Rapid Timing units |37j). and thus will be the central point in this paper. Rasnow 
in 1996 gathered these two previous results by considering a time-harmonic and spatially uniform 
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background electric field. In the presence of a splrere witli center at and radius a, conductivity 
CTi, and permittivity ei, the uniform background electric field Eq with frequency w is modified by 
adding the following field: 

f a Y {(Ji + iujEi) - {(To + iujEo) 

ho ■ X — [ -. — r — r -. r, (l-l) 

\|a;|y 2(<7i + lujei) + (aa + lujEo) 

where the index refers to the ambient medium. 

Numerical approaches have also been driven since the 70's: in 1975, Heiligenberg proposed 
a finite differences scheme to calculate the field created by the fish |52j. In 1980, Hoshimiya 
et al. use finite elements to solve this problem [3S]. The geometry of the fish is simplified by 
an ellipse and is divided into two areas: the thin skin with low conductivity and the interior of 
the body. Their aim is to optimize conductivity values to approximate as better as possible the 
experimentally measured field. The result is that the optimal conductivity is non-uniform, being 
higher in the tail region. Improvements of these models since then can be found in [T^l 1331 [351 [32] 
and references therein. However, the most promising technique is the use of the boundary element 
method performed by Assad in the 90's in his PhD thesis [10]. Indeed, the important feature is the 
electric potential on the skin (because it is the input for the fish), so a boundary element method 
(BEM) approach allows us to concentrate the equations on it. Moreover, the computation speed is 
enhanced because the number of nodes is dramatically reduced. The equation considered is here 
Aw = on the exterior of the body with Robin boundary conditions on the skin |51j : 

"-«£ = *■ <'-^) 

where ip is the potential inside the body and ^ — h{(7o/as) {h being the skin thickness and tJs 
(resp. (To) the skin (resp. water) conductivity) is the effective skin thickness. 

Let us mention that there are other kinds of simulations, based on a more empirical approach, 
determining an equivalent electric circuit [18( 119] or an equivalent multipole |22; . 

The aim of this paper is to derive a rigorous model for the electrolocation of an object around 
the fish. Two problems arise: the direct problem, i.e., the equations involved and their boundary 
conditions, and the reconstruction itself. For the direct complex conductivity problem, we show 



using layer potential techniques the validity of ( 1.2 ). We also generalize formula ( 1.1 1 to the case of 



a non-uniform background electric field, taking into account the distortion induced by the body of 
the fish, and with any shape of the target. For the inverse problem, little is known in the complex 
conductivity case [TS]. Here, we take advantage of the smallness of the targets to use the framework 
of small volume asymptotic expansions for target location and characterization [B] [7j . However, 
since the electric current is generated by only one emitter at the tail of the fish (the electric organ) 
and measured by many receptors on the skin, standard non-iterative algorithms such as MUSIC 
(standing for Multiple Signal Classification) cannot be applied for location search. In standard 
MUSIC, the data (called multistatic response matrix) form a matrix and its singular value decom- 
position leads to an efficient imaging function by projecting the Green function of the medium onto 
the significant image space [21111 [HI [13 IHl [231 [30] • Here, roughly speaking, one has only a column 
of the response matrix. However, using the fact that the electric current produced by the electric 
organ is periodically time dependent with a known fundamental frequency, we extend MUSIC ap- 
proach to multi-frequency measurements by constructing an efficient and robust multi-frequency 
MUSIC imaging function. We perform numerical simulations in order to validate both the direct 
model and the multi-frequency MUSIC algorithm. We also illustrate the robustness with respect 
to measurement noise and the sensitivity with respect to the number of frequencies, the number 
of sensors, and the distance to the target of the location search algorithm. Finally, in the case of 
disk- and ellipse-shaped targets, we provide a method to reconstruct separately the conductivity, 
the permittivity, and the size of the targets from multi-frequency measurements. We mention that 
this is possible only because of multi-frequency measurements which yield polarization tensors with 



complex conductivities. It is well-known that polarization tensors for real conductivities cannot 
separate the size from material properties of the target [7] . We also mention that the use of different 
values for the frequencies is more crucial for the material and size reconstruction procedure than 
for the location step. In fact, in the presence of measurement noise, location with N realizations 
with one frequency is comparable to the one with TV different frequency values. 

The paper is organized as follows. In section[2] the model is set up and the equations governing 
the electric field are rigorously derived. Using layer potential techniques, the boundary condition 
(1.2) is recovered. In section pi a small target is located using multi- frequency measurements. For 



this purpose, a dipolar approximation is derived before the analysis of the response matrix. Finally, 
numerical simulations are performed in section HI due to the presence of a hyper-singular operator, 
a particular attention is paid to the numerical scheme. Reconstructions of the electromagnetic 
parameters and the size of disk- and ellipse-shaped targets are also provided. 

2 The forward problem 

The aim of this section is to formulate the forward problem. After the setup of the problem in 

subsection 

subsection 



2.1 the boundary conditions are announced in subsection |2.2| before being derived in 



2.2.2[ Existence, uniqueness and a useful representation lemma for this derivation are 



proved in subsection 2.2.1 



2.1 Non-dimensionalization and problem formulation 

In this subsection, we derive the equations governing the electric field. A formal explanation of 
the electroquasistatic (or EQS) formulation is given, and the setup of the problem is then non- 
dimensionalized . 

Partial differential equations of the problem 

The electroquasistatic (or EQS) formulation is a low-frequency limit for the Maxwell system in 
three dimensions. In the frequency domain, this latter is given by 



(2.1) 



where E is the electric field, B is the magnetic induction field, p and j are the free charges and 
currents, w is the frequency, /i is the magnetic permeability, and e is the electric permittivity. 
Moreover, in a medium of conductivity a the Ohm's law connects the electric field to the induced 
current density {ji = aE) so the current density can be decomposed as: 

j = aE + js, 

where js is a source of current (in our model, it comes from the electric organ). Then, taking the 



V -sE^ 


= P, 


V B- 


= 0, 


V X E^ 


= -icoB 


/i 


= j + iu! 



divergence of the last line in (2.1 1, we have: 

V ■{a + ieuj)E = -V -js- (2.2) 

The EQS approximation consists in considering the electric field as irrotational because the mag- 
netic field variation is negligible. A sufficient condition for that is given by [15] : 

^ « 1, (2.3) 



where imax is the maximal length of the problem and Anjju the minimal wavelength. Here, we can 
take Lmax ~ Ini because the range of electrolocation does not exceed two body lengths |37_. In 
the water, the minimal wavelength is given by 

1 



[IE 



where /i « /xq, e « 80eo a-nd Wmax is the maximal frequency emitted by the fish, which is of the 
order of lOkHz. Thus, the fraction in (2.3) is of order 10^"*, so the EQS approximation is very well 
suited for our situation. 



Going back to the equation of the electric field (2.2), we can now use the fact that E is 



irrotational to state that it is derived from a potential scalar field u. This finally leads us to the 
following equation: 

V • (cr + iew)Vu = "V- js- (2.4) 

To conclude, taking into account the slow variation of the electric field leads us to consider 
a complex conductivity instead of a real valued one. However, for the rest of this section, the 
imaginary part of this conductivity will be neglected; indeed measurements on a Gnathonemus 
petersii showed that the permittivity of the skin, the body, and the water are very small compared 
to their respective conductivity |19[ I40j . Thus, this EQS approximation will be used only in the 
presence of a target: it will be detected by the phase shift induced by its complex conductivity. 

Non-dimensionalization 

We wish to perform an asymptotic analysis of the equations. The first step consists in the identi- 
fication of the different scales of the model problem. The electric potential u, the variables x and 
Lo, and the parameters a and jg can be written as follows: 



U = VnU , X — Lx , LO ~ 






where Vq is the voltage produced by an electric organ discharge (EOD), L is the length of the fish, 
wq is the fundamental frequency of the EOD, (Tq is the conductivity of the surrounding water and 
/o is the current intensity inside the electric organ. Moreover, anticipating the next subsection, 
the conductivity of the body and the skin play an important role in the shape of the electric field. 
Thus, in the list of parameters we add the conductivity of the body ah^ the thickness of the skin h 



and its surface conductivity E. The orders of magnitude of these parameters are found in Table 2.1 



Quantity 


Order of magnitude 


Reference 


Vo 


10 mV 


Ui.Sl 


L 


10 cm 


37 


Wo 


1 kHz 


.37 


cro 


100 /iS-cm-i 


.34 


h 


I mA 


.14 


CTfc 


I S-m-^ 


.40. 


S 


100 /iS-cm-^ 


19 


h 


100 ^m 


52 



Table 2.1: Orders of magnitude of the physical quantities involved. These are only scales and not 
the exact values measured in the cited references. Here S is Siemens (IS* = lA/iy). 



These n — S quantities involve r 
Buckingham-Pi theorem, we need n - 



= 4 fundamental units of the SI system, so according to the 
r ~ A nondimensional quantities. The first one can be found 



by rewriting the equation (2.4 1 in terms of the nondimensional quantities {x' , k,u',j[.] 



V,.-kV,,u' = ^V-j-;. (2.5) 

The niultiphcative term in the right-hand side of the previous equation is not important as the 
equation is hnear. The three other nondimensional quantities come from the parameters of the 
skin and the body of the fish: 

fc, := ^ ^ 102, fc, :^ ^ ^ 10-2, 5:^^^ IQ~\ 

In other words, in nondimensional units, kjj (resp. k^) is the body (resp. skin) conductivity and 5 
is the skin thickness. 

To conclude, omitting the prime symbol for the sake of clarity and denoting by fcf,/ the source 



term in equation (2.51, the governing PDE is the following 

V • kVu = kbf, (2.6) 

where k is piecewise constant, being equal to 1 in the water, kh inside the body of the fish and kg 
in the skin. These domains are going to be made precise in the next subsection. 

For the sake of simplicity, from now on, we only consider the model equations in two dimensions. 

2.2 Boundary conditions 

In this subsection, we derive the appropriate boundary conditions associated with the presence of 
a very thin and very resistive skin. Robin boundary conditions will be found after an asymptotic 
analysis of the layer potentials involved. 

The setup is as follows: the body occupies a fixed smooth open set 51(, and the skin with 
constant thickness is described as: 

ils :== {x + tv{x), X e 95^6, < t < (5}, 

where v is the outward normal unit vector. Let us also denote by ^ the effective thickness defined 
by Assad |51j : in our variables it is given by 

The source of the electric field is a sum of Dirac functions: 

m 

where, for \ < j < m^ Zj G flh and / satisfies the charge neutrality condition 

m 

J2a,=0. (2.7) 



Although condition (2.7) is the physical condition in our model, we will show how to modify the 



derivations and the results of the paper in the general case. An illustration is given in Figure [271 



Our main purpose here is to investigate the behavior of the solution of (2.6) with 



ks a X G ils, 
k{x) = { kb iixe nt, (2.8) 

1 otherwise. 



c 



= <5S 



(^b 




Figure 2.1: Setup of the problem. The conductivities are non-dimensionahzed so that ctq — 1. The 
body fib is represented in grey and the skin fl^ is represented by its bold boundary. The sources 
/ are given by the two dots. 



where kg ^ 1 and k^ ^ kg, in the following asymptotic regime: 

5 



e 



, ^ is fixed, 



0, and fcfc — )■ oo. 



In order to make this dependence clear, let us denote such a solution by ug^k^- Adding a far field 



condition (essential for uniqueness, see subsection 2.2.11, it is the solution of 



V • kWus^k^ = kbf, 



X e 



\U5M 



0{\x\ ), \x\ — > oo uniformly in x, 



(2.9) 



where x := x/ \x\ and k{x) is given by (2.8). Note that if assumption (2.71 does not hold, then the 



far field condition should be replaced with 



.^ . (A6 + l/2)(A, + l/2) 

3 = 1 



'2^(Ab-l/2)(A, -1/2) 



= 0{\x\ ), |a;| — >■ oo uniformly in x, (2.10) 



where the parameters As and Xb are given by 



A. := 



2ikg 



1 J , kg+kb 

— and Xb := 



2{kg - kb) ' 



(2.11) 



The far field condition ( 2.10 1 will be explained later. We will compute the first-order asymptotic 



wo,oo and see that it is the solution of the following system: 

A-uo^oo = /, 
A-uo,oo = 0, 

duQr 



■"0,00 |_|. ~ "0,00 

duQoo 



^ 



dv 



dv 

|"0,oo| 



X E fib, 

xeR^X^b, 
X e drib, 

X £ dflb, 



(2.12) 



0{\x\ ), |a;| -^ oo, uniformly in x. 



Note that in the limiting model (|2.12|), the role of / is to fix the potential uqoo\_ on dflb- On the 



other hand, if assumption (2.7) does not hold, then the boundary condition on g^ 
replaced with 



should be 



diy I- |au.i^"^ 



a, 



j=i 



To be more precise, we will prove the following theorem: 



Theorem 1. There exists a constant C independent of 5 and ki, such that the following inequality 
holds for S and 1/fcb small enough: 



\\US.kb 



"O.ooIIlo 



<C{S+^ 
kb 



(2.13) 



where us,kb 'J^'^ uq.oo o,re the solutions of (2.9) and (2.12), respectively. 



In a first part, we will analyze equation (2.9) and show that there exists a unique solution 



that can be represented as the sum of a harmonic function and two single-layer potentials. In a 
second part, we will perform asymptotic analysis of these layer potentials in order to show that 



the limiting function is solution of (2.12 1. This latter part is to adapt the work done by Zribi in 



his thesis [53 and by Zribi and Khclifi in [29 . 



2.2.1 Existence, uniqueness, and representation of the electric potential 



In this part, we will first prove the uniqueness of the solutions of (2.91 and then we will derive a 
representation formula, which will give us the existence of the solution. For the moment, 5 and k^ 
are fixed, but we suppose that: 

ks<l< kb. (2.14) 



Uniqueness 

The uniqueness comes from the second line of ( |2.9[ ) [7]. Indeed, let v — ui — U2, where ui and 1*2 
are two solutions of (2.9) and let us show that v = 0. From (2.14) we get, for R sufficiently large 



(so that the ball with center and radius R encompasses Sis): 



iVwl' < -^ 



|x|<_R 



k{x) |Vu|' 



|x|<fi 



dv 
ks .l\x\=R dv 



iVi^r <o. 



s J\x\>R 



Here we have used the fact that Vu G _L^(M^ \ fi^), which holds as a consequence of the far field 
condition. A unique continuation argument shows that |Vf | = in M^ and thus v is constant. 
Then, using the fact that w — > as |a;| —>■ 00, we have w = 0. 

Existence and representation 

The existence is given by a representation formula decomposing the solution into a source part 
and a refraction part. This refraction part implies layer potentials on the boundaries of the body 
and the skin. Let us define them explicitly and give some well-known results. First, let us define 
the following boundaries: 

Tb := dilb and T, := dn^ \ Tb. 

In the following, the index f3 stands for the subscript b or s. The single- and double-layer potentials 
on F/j are operators that map any (p G L'^(Tp) to Sj3(p and T^j^ip, respectively, where 



Sp :— Sr^ with 5r<,5 := / G(- — s)ip{s)ds, 
Dp :=: Vti, with Vrip := / 7;— (• - s)ip{s)ds, 



where G is the Green function for the Laplacian in M^: 



G(x):=— log|a:|, x ^ 0. 

ZTT 



(2.15) 



For ip e L^(r^), the functions SpLp and "Dpip are harmonic functions in M^ \ F^; their singularities 
hold on F^. To describe these singularities, we define, for a function w defined in M^ \ F^ and 

X e F^: 

w{x)\^ := Ihn w{x zLtiy{x)), 



dw 
dv 



(x) 



lim Vwfa; ± tvix)) ■ vix). 



Across the boundary F^, the following trace relations hold [7j: 

dSpip 



±2^ + K^} ) f, 



dv 



Vp^\^^ \^-I + ]Cfi\v 



(2.16) 



&Dpifi 



dv 



dVpip 



dv 



Here, the operator Kp and its L'^-adjoint /Co are given by 



(/C^(^)(x) := 
(/C^(^)(x) := 



1 

2^ 
1 

2^ 



p.v. 



p.v. 



{s-x) ■v{s) 

— ^2 — V{s)ds , a; € F/3, 

\x~ s\ 

{x - s) ■ v{x) 



where p.v. stands for the Cauchy principal value. From (2.161 it follows that the following jump 
formulas hold: 



dSpip 



dv 



dSpip 



di 



ip and Vpifil - Vi3ip\_ 



-if. 



The following invertibility result is useful PU 37] . 

Theorem 2. Suppose that Tp has Lipschitz regularity. Then the operator XI — 1C% is invertible on 
L'^iVp) := {ipe L^iTp) : f^^ip = 0} if \X\ > 1/2, and for X e (-oo, 1/2] U (1/2, +oo), A/ - /C^ is 
invertible on L'^(Tp). 

With these essentials tools, we can now prove the following decomposition formula in the same 
spirit as in [77] : 



Lemma 1. The solution of problem {2.9) can be written as 

u{x) = H{x) + {SsiPs){x) + {Sb(pb){x), 
where 

m 

and the pair (ips^cpi,) G L^(Fs) x _L^(Fh) is uniquely determined by 

[XJ - JC^jips J. — = ^-, X G F^, 

dv dv 

[Xbl - K.bJ'Pb ^ = -7^, X G Tb. 

dv dv 



(2.17) 
(2.18) 



(2.19) 



Here, Xb and Xs are given by (2.11). Moreover, the decomposition {2.17) of u into a source part 
H and a refraction part Ssifis + SbPb is unique. 



Proof. The system (2.9) is equivalent to the following transmission problem [T]: 



Au^f, 



X e 



u\+ - u\_ =0, 



du 
dv 
du 
1^ 



h 



du 
dv 
du 
dv 



0, 



0, 



'\(r6ur,), 
xeTfeur,, 

X e r., 



|m| = 0{\x\ ), \x\ — > cxD, uniformly in x. 



The existence of a solution {(ps,^Pb) to (2.191 comes from the fact that |As|, |Ab| e (1/2, +00) and 
Theorem [2J On the other hand, the functions Sg'^s and Sb'Pb are harmonic in ilt,, and according 
to the definition of H , we have Au = / in ^b- In ^s and M^ \ fi^ U Vtb, all these functions are 
harmonic so we have Au = 0. The trace relations on Tb and Ts are then given by the singularities 



(2.16) of 5s and Sb (see [7]) since H is smooth away from the points Zj. Finally, all the se functions 
are controlled by |a;| when |a;| — >■ 00. In this way, the existence of a solution to (2.9) is proved. 



To prove the uniqueness of the decomposition, let us take Lps and (p[, such that 

H + Ss^s + Sbipb = H + Ss(ps' + Sbip'b- 

Then, Ss{^s~'^s ) = Sb{iPb — ^b) is harmonic in figUilb, which gives by the jump formula (pb = f'b- 
Finally, applying once more the jump formula, we have ips — (ps ■ □ 



We now check the far field condition stated in (2.10). Recall that /C(,(l) = /Cs(l) = 1/2- From 



dSbf b 



= / </56, 



dSgifs 

dv 



= 0, 



dH 

dv 



dH 
dv 



= Z!"j' 



by taking the average of the two equations in (2.191 on Fg and F;,, respectively, we find that 



J 'P'^ == (5Z "j 



(Afc + 1/2) 
(A,-l/2)(Ab-l/2) 



and 



r/''-A,-l/2' 



and therefore, from the representation formula ( |2.17[ ) it follows that 

(A5 + l/2)(A, + l/2) 



WAfc, 



(E' 



'27r(A5-l/2)(A, -1/2) 



log l^l 



0{\x\ ), |x| — )■ 00 uniformly in x. 



Note that in the limit (5^0 and kb 
— 1/2, As — > —1/2, and therefore, 

l'"o,oo| — 0{\x\~ ), \x\ — > 00 uniformly in x. 
For the system (2.121, Lemma fl] yields the following result. 



00, the far field condition above and (2.111 yield A;, -^ 

(2.20) 



Lemma 2. Assume that ^2.1\ ) holds. The solution of problem ^2.1^ can be written as 

u{x) = H{x) - -{Sb(p){x) + {T)bip){x), 



(2.21) 



where H is given by (2.18) and ip G L^lTb) '■= {4> G ^^i^b) : /p = 0} is given by the following 
integral equation: 



9^ 



dv ' 



The decomposition (2.21) of u into a source part and a refraction part is unique. 



In the general case, (2.22) should be replaced with 
1 /U ...\ dVbip 



e V2 



-I-ICt 



f- 



dv 



dv 



fjE"^' 



X eVb 



Note that since /p ^ = Ylj o^i' the far field condition (2.20) is satisfied in the general case. 



The proof of this lemma involves exactly the same arguments as in the previous one: jump 
formulas applied to the operators. 



The decomposition formulas (2.17) and (2.21) will be essential in the next part to show that, 



at the first-order, us^k^ converges to uo,oo- 

2.2.2 Asymptotic expansion of the electric potential for highly resistive skin and 
highly conductive body 

In this part, we will use the decomposition formula for us^k,, and compute asymptotic expansions 
of the refraction part. The limiting solution will then be uo,oo- This latter is well defined if the 
limits 5 — )• and fc;, — >■ oo are independent, so we must seek the two following limits: 

lim lim u^ fcj_ and lim lim ug^^^ 

fcfe— »oo (5— >-0 ' {5— vO A,'b— foo 

and show that they are the same. Zribi |531 chapter 3] studied the case when fcb remains fixed, 
with non-uniform thickness of the skin $7^; the limit uq^i is the solution of the system: 



Auo,i = /, 
Auo,i - 0, 



^o,i| 



uo.il 



duo; 



dv 



duo,i 



di 



= 0, 



duQ- 



x£R^\ Ub, 

X e dfib, 
X e fib, 



(2.23) 



|uo,i| = 0(|a;| ), |a;| — > oo, uniformly in x. 
Here, we will follow the same outline for the proof: first we will remind the asymptotic expansions 



of the operators involved in (2.19), and then we will match the asymptotic expansions for (fs and 
Asymptotic expansions of the operators 



In the decomposition formula (2.17), H is independent of S and fc^; we just have to analyze the 



dependence of (ps and ifb- Remark that from (2.191 

• the dependence on kb is carried only by Xb since Sb and IC^ depend only on the shape of fib] 

• the dependence on 6 is carried by A^, Sg, /C* and d/dv{x) for a; £ F^. 

In this subsection, we will focus on the asymptotic expansions of the operators (the limits of As 
and Af, are obvious) . They have been performed in |21 153| ; in order to apply this proof, we first 
need some assumptions. 

Suppose Tb is defined in the following way: 

Tb:=9{dB), 

where 5 is a C^'"^ diffcomorphism of the unit sphere dB := dB{0, 1) for some 77 > 0. Moreover, we 
suppose that the function Xg : [0, 2tt] -^ R'^ defined by 



Xg^g 



cost 
sinf 



10 



is such that |X'(i)| = 1 for all t G [0, 27r]. Thus, Xg is a C^^'' arclength counterclockwise 
parametrization of Tf,. Then the outward unit normal to Qi,, vix) at 2; = Xg(t), is given by 

where R-2l is the rotation by —-k 12. The tangential vector T(x) at x = Xg{t) is defined by 

T{x)^X'g{t), 

and X'g{t)±Xg{t). The curvature t{x) at a; = Xg{t) is defined by 

X'^it)^TixMx). 

Let ^5 be the diffeomorphism from r^ onto Tg given by 

^six)^x + Si^{x). (2.24) 

With these assumptions, the following regularity result holds [31] : 
Theorem 3. Let 77 > 0. Let, for a Lipschitz function g e C°'^ {dB, M^), 

9{x) - g{y) 



Idsig] ■■= inf 

x^yedB 



x-y 



Introduce the set Aob of admissible diffeomorphisms of the unit sphere: 

Aqb ■■= {g e C\dB,R^),laB[g\ > 0} . 

Then, for any integer m > 0, the operators S and D defined on {C"'''''{dB,R'^) D Aqb) xC™"i'''(aS) 
((^C"'-v(^dB,R'^)nAdB) xC™^''(9B), respectively) toC"''''{dB) by 

S[g,(p]{x) ■.^Sg(OB){f ° g^^) ° g{x),x e dB, 

D[g,(p]{x) ■.^Vg(ijB){ip o g-^) o g{x),x G dB, 

are jointly analytic with respect to their variables g and ip. 

Moreover, we have explicit formulas for the derivatives with respect to the variable g |31| . 
Then, we have the following asymptotic expansions [HI [53] : 

Proposition 1. Let ip G C^''^{Ti,) and ip G C^'^'(Ts) for some r] > 0. Then, we have the following 
asymptotic expansions for x G Fh ; 

(/C:vi) o 'i'six) - ICli^(x) + 5K.\^^i,{x) + 0{6^), 
^ ° ^^six) - Q/ + JCl^ v{x) + STZM^) + 0(<5i+"), (2.25) 



dv 



{x) ={-\l + ^l) ^(^) + 5Ciij{x) + 0{5^+'^), 



where ijj := ijj o ^!g, ^ig being defined by (2.24-), ^.^d 



ici^^ipix) = t{x)1CIiIj{x) - ICI{tiI;){x) ^j^j^i^) + ~^(^) 

TZbipix) = t{x) {2^ + 1^1] "fii^) dj^^^''' 

C^ix) = Q/ - ICl^ {r^){x) + ^(x), 
where d/dt is the tangential derivative in the direction ofT{x) = X' o X'7^{x). 



(2.26) 
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Note that, according to Theorem pi the constants in the 0((5^+'') terms depend on ll^H^a,^. 
Moreover, since the thickness of fig is uniform, we have vo^ig^x) = iy{x), and a Tayfor expansion 
of H gives, for x G Ff,: 



|_ o ^s{x) ^^ix)+ ^^(^) • [D'H{xHx)] + 0{6^), 



(2.27) 



where D^H denotes the Hessian of H. 



Asymptotic expansions on the layers 

In order to prove Theorem [ll we wih first show the convergence on the layers (see next lemma). 
Then, in the next subsection, we will extend the domain of validity by application of the maximum 
principle. 

The following lemma holds. 

Lemma 3. There exist constants C and C independent of 5 and fcf, such that the following in- 
equalities hold for S and 1/fcf, small enough: 



\\u5,k, -uo,oo|li=o(r,) <C (6+ —j , 
\\u5,k, -uo,oo|lioo(r,) <C" (6+ —j , 



(2.28) 



where us^k^ and uq.oo are solutions of (2.9) and {2.12), respectively. 



Proof. Only the first limit will be shown, the second one being very similar. For this purpose, we 
must show that the limits (5 — > and fc^ — >■ oo are independent, i.e., they commute. First, let us 
compute the limit of us^k^ when (5 — > 0, and then the limit fcf, — > oo (which will be much easier). 
Then, we will invert this process. 

This first limit is the main problem in ^53t chapter 3], except that, in that study, fc;, = 1 and the 
thickness of fi^ is non-uniform. According to theorem pi the formulas in ^31j, and by composition 
with the regular diffeomorphism ^5 from Tg to Fh, we have 



Ss(fs — Sb(ps 



I + ICb] ips - Sbinps) 



<C5\ 



C2.'7(r6) 



where Lps '■— ifs o "^s- Hence, with the help of the decomposition formula (2.171, we have the 
following asymptotic expansion uniformly on F^: 



usMi^) = H{x) + Sbiipb + ^s)ix) +S 



--I + JCb ] ipsix) ~Sb{TLPs) 



{x)+0{5^), (2.29) 



We now look for expansions of the functions (ps and Lpb when (5 — ;• that will be re-injected in this 
equation. Using Proposition II| and (2.19), these functions are solutions of the following system: 



-'fs 



I + JC*b] {ips + 'Pb) + 5 



^(1), 



OH 
dv 



+ 0(<5i+'') = 
5[vD''Hy\+0[5''), (2.30) 



dH 
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Let us define the formal asymptotic expansions: 



^s = l^i~''> 



,(0) 



ifl"' + 6ip\ 



(1) 



fb 



-*- (-1) I (0) , X (1) 



Aiming to have the 0-order term in the expansion (2.291, here we seek for the terms of order — 1 



and 0. By substitution into (2.30) and identification of the leading-order terms in the first line, we 
get: 

/I \ / , .. , .^\ 

0, 



(i/+>c:)(ri- 



'-^i-' 



so that, by Theorem [2] we have: 



v:'i-^)+^r'=o. 



(2.31) 



Let us now look at the 0-order terms; summing the two lines, we get: 



= 0, 



which gives, with the help of (2.26) and (2.3l|), 



1 


.\ 1 1 




-1] -+T 


kh 


M J 



f. 



(-1) 



(2.32) 



J-i) 



This quantity is what we need in (2.291; thus, only tps re mains to be found. This can be done 

'(1 

^b 



by identification of the 0-order terms in the first line of (|2.30|) and using the definitions of ICi and 



Kb given by ( |2.26[ ): 
1 



.(-1) 



1 / 1 



^s - + -^^-^-l][-,I + ICl]^i-'^ 



dVbip 



(-1) 



dp 



dH 
dv 



Finally, the expansion (2.291 yields 



us,kb{x) = H{x) + 



1/1 



^\k 



- 1 hS, 



-^i + ic. 



^^r^\x) + o{6). 



(2.33) 



(2.34) 



This leading-order term (denoted uo^kb) verifies (2.23) according to (2.331 and jump formulas (2.16). 



The asymptotic /cf, — > oo does not add further difficulty. Indeed, let us define the following 
asymptotic: 



By substitution into equation (2.33) and identification of the leading-order terms, we get: 



-(- 



i-icU ^i-'-'^ 



dVbif 



(-1,0) 



dv 



dH 
dv ' 



and then the expansion (2.34 1 becomes 



us. 



kAx) - H{x) - ^5,^1.-1'°) (a;) + {-ll + ICb^ vi-'^"\x)+0{S) 



(2.35) 



which is ( |2.21 ) applied on Tb according to the jump formula of Vb (2.16). Hence, according to 
lemmal2j the first-order asymptotic of ug^kt is uq^o- 
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Let us now show that the hmits 5 — )■ and fcf, — > oo commute: unhke in the previous subsection, 
we will first perform the limit fcf, — >■ oo and then the limit 5 — )■ 0. Given the fact that 



1 



Xb = -^ + - 



1 



the definition of (ps a-nd iph in (2.19) will be affected only in the second line. Indeed, with the 
following expansions: 

kb 

(0) ,1 (1) , 

Kb 

this second line becomes, at the leading order: 



dv ' 



where ips .= ifs ° "^s- With the expansion: 



(0) 



V5, 



(0) 

fb 



r'P. 



(0,-1) +^(0,0) ^^^(0.1) + 



1 



:^b 



(0,-1) 



n 



(0,0) 



Sip 



(0,1) 



the identifications (2.31), (2.32 1, and (2.33) respectively become: 



(2.36) 



<P 



(0,0) 



(0,0) 

n = 



1 /I 



CV2 



-I -JCl] p^ 



.(0,-1) 



1 
dT>b(p 



(0,-1) 



(0,-1) 



dv 



dH 



(2.37) 



(2.38) 



same expansion (2.35) 



Finally, recalling that the expansion of us.kt in (2.29) is conductivity-independent, we obtain the 



D 



Proof of Theorem [T] 



With the estimates (2.28) on the layers Tb and F^, we are now ready to prove the estimate (2.131 



on the whole space applying the maximum principle. 

For the sets fib and fig, it is straightforward: the function us^k^ ~ '^i),oa is harmonic in these 
bounded domains, so the maximum is reached on the boundaries [45 . And, since this maximum 
is dominated by (5 and l/fc;,, we have: 



lW<5,fct 



Wo, 



°o|lL~(06Un^) 



<cu 



1 



For the exterior domain, we cannot apply directly the maximum principle since this domain is 



unbounded. However, the conditions at infinity in the systems (2.9) and (2.12) allow us to have a 
similar control. Indeed, this condition tells us that 



[""^.fc, 



"0, 



oollL°°(B(0,fl)) 



= 0(i?-i) 



(2.39) 
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We take: 



^ •— 2 11^''^'='' "o^°ollL°°(nbunj ' 



and choose Rq such that, for R > Ro, the right-hand side of (2.39) is bounded by e. Then, we 
have: 



\\us,kh - WO,oo|lj 



J2\B(0,i?,o)) 



<£. 



Now, only the bounded domain B{0, Rq) \ [fli, U fls) remains, where we can apply the maximum 
principle. Thus, Theorem IT] is proved. 



2.3 Final formulation and notation 

In the previous subsections, we have performed a multi-scale analysis of the problem to identify 
the effective equations with boundary conditions. In order to make things clear, let us summarize 
the results and simplify the notation. 

The electric potential emitted by the fish is the solution of the complex-conductivity equation 
(2.6) with boundary conditions given by the system (2.12). It is easy to see that in the case of 



an inhomogeneity outside the body, these boundary conditions will not be changed because the 
asymptotics are done with the layer potentials of the domains defining the fish. 

Hence, we conclude this section by summing up the results: omitting all the subscripts, the 
electric potential u is the solution of the system 

xeR^\n, 



n 


^u = /, 


V • (1-1- (fc - 1 + ieuj)xD)^u = 0, 


^ du 


= 0, 

+ 


du 


= 0, 




u\ = 



(2.40) 



xGT, 

x&T, 

), |x| — )■ oo, uniformly in x, 

where xd is the characteristic function of the target D, k + ieoj is the conductivity inside D, uj 
is the frequency, and k and e are positive constants. Here, we have assumed that ^ aj — 0. In 
the case where it is not, we should replace the boundary condition §^| = on F with §^|_ = 
(l/|r|)^ Q!j. From now on, we restrict ourselves to the case J^j'^j — 0- Note that taking two 
points zi and 226^ close enough and ai — —oii ^ yields an approximation of a dipolc at 
(z\ -I- Z2)/2 of moment |ai| and direction orthogonal to (z\ — 22). 



3 Detection algorithm for mult i- frequency measurements 

In this section, we develop an algorithm to recover (from a single measurement) the location 
of a small object located far away from the fish. This algorithm is based on multi-frequency 



measurements, as it is explained in subsection 3.1 In subsection 3.2 asymptotic expansions will 
be carried out for the electric field in the presence of a small and distant target. Finally, the 
algorithm will be explained in detail in subsection |3.3[ 



3.1 Multi- frequency measurements 

Let us suppose that the electric current produced by the electric organ, (i.e., the source term / in 
equation (2.40l) is periodically time-dependent with separation of variables, that is 



]{x,t)^j{x)h(t), 
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where / is a sum of Dirac functions and h{t) is periodic with fundamental frequency loq. Hence, 
we set 



N 



h{t) = ^ /i„e 



inojQt 



(3.1) 



where N is the upper bound that ensures the low-frequency regime. According to the previous 
section, the electric potential u is then given by 



N 



u{x,t) ^ y^ hnU^e 



inujQt 



(3.2) 



where m„, for n — 1, . . . , N, is solution of the following system 

Au„ = /, 
V • (1 + (fc - 1 + ienLOo)xD)Vun ^ 0, 

du 



X £ 



^ 



dv 

du„ 
dv 



0, 



0, 



a; G 17, 

iR2\n, 

xer. 



xer. 



(3.3) 



|m„| — 0{\x\ ), |a;| — > oo, uniformly in x, 



3.2 A dipolar expansion in the presence of a target 

In this subsection, we derive useful formulas in order to simplify the data. For the sake of simplicity, 
and for numerical reasons that will be given in section |4j only one target D will be considered. 

The electroreceptors of the fish measure the electric current at the surface of the skin [S7] . 
Hence, from a single measurement, we can construct the Space-Frequency Response (SFR) matrix 
A, whose terms are given by 



Aln = 



dUn 

du 



dU_ 
dv 



[xi), for 1 < 71 < TV and 1 < Z < L, 



where (a;z)i<;<i are points on the boundary F and U is the static background solution, i.e., the 
electric potential without any target which does not depend on n. It is the solution of (3.3) with 
a constant conductivity equal to 1 outside the body f2. 

The first formula, given in Proposition [21 is often called a dipolar expansion; indeed, in the 
presence of a small inhomogeneity, the perturbation of the electric potential looks like the electric 
potential of a dipole O [2^ . More precisely, using exactly the same arguments as in [6l Chapter 4] 
and in [5| we have the following result. 

Proposition 2. If D := z + aB with dist(z,F) 3> 1. a <C 1 and B is an open set, then we have 

^ OGr 



Hn 



-a^VUizfM{kn,B)V, 



dV:r. 



{xi,z), 



(3.4) 



where T denotes the transpose, kn = k + ieuj^n is the (complex) conductivity of the target at the 
frequency nujo, M{kmB) = {Map{kn, B))a^p=i,2 is the first- order polarization tensor associated to 
B with conductivity kn J^: 



Mapikn,B) := / (XrJ-lCg) ^{va)vtjds{y), \r. 

JOB 



kn + l 
2{kn - 1) ' 



a,/3 = l,2, 
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and Gj^ is the Green function associated to Robin boundary conditions, which is defined for z € 

R^\Tl by 

-A,Gr{x,z) = S,{x), xeR''\n, 



Gf 



C 



dG 



B. 



1 



= 0, 



xeT, 



(3.5) 



Proof. Let 



We have 



and on the other hand, 



zvr 



H„ 



= 0{\x\ ), |a;| — )■ oo, uniformly in x. 



,du„ 



-M^U)+Vr{u„U). 



Un-U = -(fc„ - 1) / VUn ■ VGr, 
J D 

Un - Hn = -{kn - 1) / Vm„ • VG. 

Jd 



From the transmission condition 



it follows that 



Since 



Un -U - 



dUn 



= on dD, 






(3.6) 



lViI„-V[/|Uoop) KCa" 



for some constant C, provided that dist(£', 9J7) ^ a, a scaling of the integral in (3.6 1 together with 
a Taylor expansion of Gr gives the desired asymptotic expansion. Note that the approximation in 
(3.4) is uniform in I and fcn [Zj- D 



Now, we will carry on a second formula in order to simplify this equation. Indeed, the Green 
function associated to Robin boundary conditions is tedious to compute. Instead, we will post- 
process the data thanks to the following lemma which generalizes Lemma 2.15 in [6]. 



Lemma 4. Let G denote the Green function in the free space defined by (2.15). For z G M^ \ fJ 
and X e F, let Gz{x) = G{x — z) and Gr^z{x) = Gr{x — z). Then 



1, ^. ,dVr 



dG 



R,z 



dv^ 



(x) 



dG, 

dvx 



(x). 



Proof. Employing the same argument as in Lemma [2] yields 
Gr,z — —G 



^-{Srv) + Vr^ - 5r(^) - C^^rl^), 



~dG 



where (p = C g^" L-- Therefore, taking the normal derivative of the above identity and using the 



trace relations (2.161 give the result 

Hence, after a calculation of ^ 
given in Lemma [4] The modified matrix will still be denoted A 



D 



Hence, after a calculation of ^j^j , ~ ^| 1 on F, we will apply the post-processing operator 
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To conclude, the location of the target D is going to be recovered from the knowledge of the 
following data 



^"-^I'-^i-i"^ 




dU_ 
dv 



[xi), l<l<L,l<n<N, (3.7) 



which is approximately equal to 

Ai,,~a^\7UizfM{K,B)V, 



ixi,z), 



(3.8) 



when the characteristic size of the target a is small. It is worth mentioning that the polarization 
tensor M(fc„,i?) is symmetric (but not Hermitian) [7]. 

3.3 A location search algorithm 

Scholz described in \A2\ a way to recover the location of a target from multi-frequency measure- 
ments. The paper focuses on an application in electrical impedance tomography (EIT) for breast 
cancer detection; the algorithm was called "Space-Frequency MUSIC" . Indeed, it is based on the 
so-called MUSIC algorithm, which is a standard tool in signal theory for the identification of sev- 
eral signals with an additive noise [H] [T^ . It has then been applied to identify small conductivity 
inhomogeneities in [31[51[T7]. In this section, we apply a similar approach for our model. 

As we can see in formula (3.81, the rows of the SFR matrix are - to leading-order - linear 



combinations of the derivatives of dG/dvx- Moreover, one has to distinguish whether the target 
is a disk or not. Indeed, in dimension 2 and in the case of an ellipse whose semi-axes are on the 
Xi-axis and of length a and &, the polarization tensor M{k, B), for k G C, takes the form [55] 

a+b 
M{k, B)^{k 



l)\B\ 



hkb 









b+ka 



Hence, the polarization tensor is proportional to the identity matrix if and only if a ^ b, i.e., B is 
a disk; this result remains true in dimension 3 [7 . This changes dramatically the range of A: if B 
is a disk, the response matrix has rank 1 and if it is an ellipse, it has rank 2. 

For the sake of simplicity, let us suppose that B is the unit disk. The identification process will be 
based on the following fact 



Lemma 5. The following map 



A 



\n^L^{r) 



z^VC/(z)^V,|^(.,z), 



is one-to-one. 

Proof. Suppose that z and z' are points on M^ \ 51 such that A(z) 
two following functions 

x^VU{zYV^G{x,z), 

v^i : M?\n\j{z'} ^M 

xh^ Vt/(z')^V^'G'(x,z')- 
Thus, these two functions both solve the following boundary value problem 

Aw = 0, xGR^\n\J{z}\J{z'}, 

dv 

^- = </?, X eT, 

ov 

w — ^ |a;| — >■ oo, uniformly in x. 



A(z') :— Lp. Let us define the 
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Hence, by the uniqueness of the solution for this problem, we have 

VU{z) ■ \7,G{x, z) ^ VU{z') ■ Vz'G{x, z'), for all a; e M^ \ H U {z} U {z'}. 
Relying on the singularity of G(-, z) at the point z, this is only possible if z = z' . 



n 



However, we do not have access to the complete function (because there is only a finite number 



of electroreceptors on the body), and the formula for A is only an approximation, based on (3.8) 



The location of the target will then be approximated as follows. In the following we suppose for 
the sake of simplicity that Xi, . . . ,xl are equi-distributed on T. 



Proposition 3 (Space-Frequency MUSIC). Define the vector 

~g{z) := {vU{z) ■ V, (^) (a;i, z), . . . , VC/(z) • V, 

and its normalized version g = g/\g\- Then, in the limit L — > 
imaging functional will have a large peak at z: 



dG 
dV'x 

oo and a 



{xL.z)] , (3.9) 

0, the following 



I{Zs) ■■ = 



1 



\{l-P)g{z,)V 

where P is the orthogonal projection onto the first singular vector of the SFR matrix A. 



(3.10) 



Proof. First of all, let us rewrite (just for this proof) the projection P^ and the illumination 
vector g^, in order to take into account the dependence with respect to L. When L goes to 
infinity, quadrature formulas show us that 



|(/-Pi^)g^(z,)|R. -> 



iI~Pa] 



HZs) 



\Hzs)\mr) 



LHr) 



Here, Pa is the projection onto the first singular vector of the operator A^, acting on the space of 
functions that have the form (3.1 1 



'-'>-^-^^ 



du 



dU 
ft7 



where u is given by (3.2| and U is the background solution (i.e., the solution of (3.3) with xd = 0). 
In the limit a — > 0, A^ is approximated by the operator A : ft, i— > A(z)h, which is obviously of 
rank one. By theory of perturbation [25], one has therefore 



U-^a) 



A{zs 



\Mzs)\mr) 



L^r) 



(I-P) 



A(z.) 



\Mzs)\ 



LHT) 



, a — >■ 0, 



L2(r) 



where P is the projector onto the first significant singular vector of A. Then, from Lemma [5) this 
functional is zero if and only if Zs = z. 

D 

Moreover, in order to have a general algorithm which is robust with respect to the background 
solution, we will plot the following imaging functional: 



I{Zs 



1 



1 



\{I-P)g^{z,)y\{I-P)g^{z,)\ 



(3.11) 
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where g^ is defined in Proposition k3^ and g^{zs) is the norniahzation of the following vector 



r{z) 



V.(g)(x„z),...,V. 



du^. 



ixL,z) 



Numerical results will be given in section [4] 

Let us highlight the fact that in the case of a general shape B, we do not know theoretically 
what happens. Indeed, when k is real, M{k, B) is equivalent to the polarization tensor of an ellipse 
[T7] , but this is not true when fc G C because the proof relies on the spectral theorem. Here, we 
still have symmetry [7], but it is not sure if M{k, B) is diagonalizable or not. However, we will see 
in the numerical subsection [4] that the algorithm works with shapes other than ellipses and disks. 

4 Numerical simulations 

In this section, numerical results are presented in order to illustrate the multi-frequency location 
search algorithm introduced in the previous section. In the first subsection, we explain the method 
used to compute the electric field; this will be the input of our location search algorithm that will 
be performed in the second subsection. 

4.1 Direct problem 

This section is devoted to the computation of the electric field around the fish. 

4.1.1 The case without target 

The electric field U generated by the fish is the function uo,oo treated in section [2| Let us recall 
that it is the solution of the following system: 



U\ 



U\ 



A 


u^f, 








xen 


AC/ = 0, 








xeR^\n 


dU 


= 0, 

+ 








xeT 


dU 


= 0, 








xeT 


|( 


7\^0i\x\ 


-'), 


\x\ -> 


CX), 


uniformly in x 



(4.1) 



Numerical simulations will be done using a boundary element method (BEM). Indeed, we need 
accuracy on the skin of the fish, and the jumps at the boundaries are too difficult to handle with 
a finite element method. Moreover, it reduces the number of discretization points, resulting in a 
much faster algorithm. 

This BEM simulation relies on the representation formula for U in terms of the layer potentials. 
From Lemma[2j we have U = H+Srip+'Dr'P, where AH = / in the whole space, and the potentials 
are solutions of the system: 



(4.2) 



Note that we have changed a little bit the notation, in order to be able to test the case ^ = 0. On 
smooth domains, the operator /Cp is easy to handle because its kernel has integrable singularity, 
whereas the operator d'Dr /dv is an hypersingular operator. Thus, one has to perform a integration 
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by parts in order to regularize it: for two smooth functions Vi and V2, we have (for example from 
[?51 Theorem 1] and |H Theorem 6.15]): 



dv 



■ V2 



G{x — y)curlrwi(x) • curlrW2(y) ds{x) ds{y), 



(4.3) 



r ^T 



where curlr is the surface rotational, defined in the following way in dimension 2. First, let us 

define the vector: 

/ di\ 

curlpi; = „_ 

ov 

where v is an extension of i; into a neighborhood of F, i.e., v{x) — v (V{x)) with the local projection 
V onto F. Then curlp is defined by 

curlru(a;) :— v{x) ■ curlp-D(a:). 



In our context, this can be made much easier. Recalling the notation of subsection 2.2.2 
have 

Thus we have, for a; G F, 



we 



F 



X = X(t) 



curlrw(x) = vi{x)- — (x) - iy2ix)^ — (x) 
0x2 oxi 

^ x^(t)^(x(t)) + x;(t)^(x(t)), t = x-\x), 



Hence, denoting by v' the curvilinear derivative of v on F, formula (4.3) becomes 

f dVyv, f , , 



This enables us to derive a BEM formulation of the system (4.2); however one has to perform it 



with Pi elements instead of simple Pq elements in the case of ^ = 0. 

The discretization process is classical [43 . We only precise that the equation is penalized in 
order to handle the condition at infinity, which fixes an additive constant. To conclude, let us 
mention that this boundary element formulation can be extended to the three-dimensional case 
(see PI). 
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4.1.2 The case with a target 



In this subsection, we derive the modification induced on t he sy stem (4.2 1 in the presence of a 
target Z? d M^ \ fi of (complex) conductivity k. The system (4.1 ) becomes: 



du 



t 


\u — f 


Au = 


du 


= 

+ 


du 


= 


u — u _ — 


, du 
+ 9^ 


= 



X e 



l^\{U^dD), 
xer, 

xer, 

x£ dD, 
xe dD, 



(4.4) 



|u| = 0{\x\ ), \x\ — >■ oo, uniformly in x. 



Thus, u can be written as 



u{x) — H[x) + Srip{x) + VTifix) + Sqd4>{x). 
The absence of Vqd is justified by the continuity across the boundary of D. From the jump 



formulas (2.16), the conditions on the boundaries F and dD given in (4.4| leads us to the following 
system: 



i-- 



r 



<'^)y^-liSon<P)\r 



dH 



du 
d d dJ-f 



xeF, 
xe dD, 



(4.5) 



dD 



where 



System (4.5 1 can be rewritten as follows: 



r 



fc + i 

2(fc-l) 



dH 
dv 
dH 



dD' 



with 



/ 






( dSr 


V 


\ dv 



/cf + c 



e 



drPy 
~d7 
dVr 



dSoD 



dD 



dv 



dD 



dv 



{XI^ICId) 



The BEM formulation is then also classical, because the only difficulty is due to the hypersingular 
operator in the upper left term. Hence, we discretize ip G Li^lT) with Pi elements and G L'^{dD) 
with Pq elements. 
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(a) 



(b) 



Figure 4.1: Isopotentials of the case described: (a) global overview (b) zoom on the target. 



4.1.3 Direct simulations 

In this subsection, we present some numerical simulations of the direct problem. We approximate 
the shape of the fish by an ellipse with semi-axes of lengths 1 and 0.3; the electric organ is a dipole 
in the a;i-direction of moment 1, placed at zq = (0.7,0) and the impedance is ^ = 0.1. A ball of 
infinite conductivity (more precisely, with a — 10^° and e = 0) and radius r — 0.05 is located at 
(1.5cos(7r/3), 1.5sin(7r/3)). Figure 4.1 shows the isopotentials. In Figure 4.1 (b) it can be seen 



that the isopotentials avoid the target since it is of infinite conductivity. 



4.2 Target location 



In this subsection, we show numerical target location results using the imaging function (3.11). 



With the same parameters used for Figure |4.1| for the fish, but with a small target of electric 
parameters a = 2 and £ = 1, we obtain the imaging functional plotted in Figure 4.2 (a). We use 
10 frequencies equidistributed from 1 to 10. In Figure 4.2 (b) and (c), we have tested other shapes 
for the target. 



Stability estimates with respect to measurement noise 

Let us first notice that, in the absence of noise, the number of used frequencies does not change 
significantly the image. Indeed, we can see in Figure |4.3| that we can recover the location of the 
target with only one frequency. 

Let us now consider the effect of measurement noise on the performance of the location search 
algorithm. We add to the entries of the matrix A defined in (3.8 1 independent Gaussian random 
variables of mean and standard deviation 



Vc 





I du„ 




dU 


^ 


ix,) 


max 




- 


- 




;,n 


\ dv 


+ 


dv 


+) 





The parameter C, is the relative strength of the noise, and will be given in %. Figure 4.4 shows 
that increasing the number of frequencies stabilizes the image. 

More quantitatively, we have computed the empirical root mean square location error (between 
the exact location of the target and the maximum of the imaging functional), for N^ = 250 trials. 
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Figure 4.2: Detection (left) of the target with the SF-MUSIC algorithm, for different target shapes 
(right). Here, the number of used frequencies is 10, equidistributed from 1 to 10, and there are 64 
equidistant sensors on the fish. 



Here, the same target as in Figure |4.3| is considered. Results are shown in Figure |4.5[ 

A natural question is whether taking different values for the frequencies plays a role. In Figure 
4.6 we use the data obtained by 100 trials for 1% of noise, 64 sensors, and a single frequency 



equal to 1. Figure 4.6 shows that the values of the frequencies do not play a crucial role in the 
location procedure. In fact, the location result is similar to the one in Figure [44} However, from a 
practical point of view, using simultaneously N different frequencies yields a faster robust location 
procedure than repeating N times the data acquisition procedure with the same frequency. In 
subsection |4.3[ we also identify the more fundamental role of the values of the frequencies in the 
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Figure 4.3: Target detection in the absence of noise, with only one frequency equal to 1. Here, the 



target is a disk with center (1.5cos(7r/3), 1.5sin(7r/3)) and radius 0.05, like in Figure 4.2 |a); the 
number of sensors is the same. 




0.3 0.4 0.5 0.6 0.7 




Figure 4.4: Influence of the number of used frequencies on the stability. Here, the same target as 
in Figure 4.3 is imaged with 1% of noise and 1 frequency (left), 100 frequencies equidistributed 
from 1 to 100 (right), with 64 sensors. The disks plot the exact position, and the squares plot the 
location of the maximum of the imaging functional. 



characterization procedure. 

The number of sensors is also crucial in the stability of the algorithm. Figure |4.7| compares the 
root mean square location error with 100 frequencies equidistributed from 1 to 100 for 64 and 8 
sensors for different measurement noise levels. 

The same type of statistics is possible for the detection as function of the distance between the 



fish from the target. In Figure 4.8 we have plotted the root mean square location errors, with 
15 frequencies equidistributed from 1 to 15 and 5% of noise, for disks with radius 0.05 placed at 
(tcos(7r/3),tsin(7r/3)) for t = 1,1.5,2,2.5, and 3. 



4.3 Target characterization 



Once the target is located, one can use (3.8) to estimate the electromagnetic parameters and the 
size of the target. Assume that the target is a disk of radius a, placed at z. From (3.8) it follows 
that a^{kn — l)/(fcn + 1) can be estimated for 1 < n < iV from the measurement matrix A. Here, 
kn — k + ieojQn with wq being known. Let t^ be the estimated values of a^{kn — l)/{kn + 1) 
from A. To characterize the target and approximate its size, one minimizes the following quadratic 
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Figure 4.5: Influence of the number of frequencies on the root mean square location error for 250 
trials. Here, the horizontal axis is for the measurement noise level in % and the vertical axis is for 
the root mean square location error. 
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Figure 4.6: Influence of the values of used frequencies on the stability. Here, the same target as 



in Figure 4.3 is imaged using the data obtained by 100 trials with 1% of noise, 64 sensors, and 
frequency equal to 1. The disks plot the exact position, and the squares plot the location of the 
maximum of the imaging functional. 
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Figure 4.7: Influence of tlie number of sensors on the root mean square location error for 250 trials. 
Here, the horizontal axis is for the noise level in % and the vertical axis is for the root mean square 
location error. 
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Figure 4.8: Influence of the distance to the fish on the mean square location error for 250 trials. 
Here, the horizontal axis is for the distance to the fish and the vertical axis is for the root mean 
square location error. 
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misfit functional: 



E 

l<n<7V 



a2(A:„-l) 



„est 



1 



(4.6) 



over fc,e, and a. 

Table [4J] gives the result of the optimization algorithm for a disk-shaped target with center cen- 



ter (1.5cos(7r/3), 1.5sin(7r/3)) and radius a*'"°. The electromagnetic parameters are (o-*''"°,e''"°). 
The initial guess is a'"'* = 0.01, a™'* — !,£'"'* = 1. The data is collected for 100 frequencies 
equidistributed from 1 to 100. The reconstructed results are accurate. 
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2.0017 


0.06 
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0.04 
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2 


0.0404 


2.9614 


1.9806 



Table 4.1: Target characterization by minimizing the quadratic misfit functional (4.6) using data 



collected for 100 frequencies equidistributed from 1 to 100. Here, true: true values, est: estimated 
values. The initial values are a'"'* = 0.01, cr'"" = 1,^""' = 1. 



When the target is an ellipse, the measurement matrix A may not be sufficient to characterize 
the electromagnetic parameters and the size of the target. At least two different positions of 
the fish (or equivalently two different locations of the target in the fish frame of reference) are 
needed in order to generate non-parallel dipole directions VC//|VC/| at the location z of the target 
and consequently lead to the extraction of the polarization tensor M(fc„, D) of the ellipse-shaped 
target D. Consider two target locations zi and Z2 m the fish frame of reference. Multi-frequency 
measurements lead to two SFR matrices, Ai^ and Ai^ with 1 < I < L and 1 < rt < A^. Define 
the following linear application from the set Jv[ of complex symmetric 2x2 matrices to C^^ 



/ 



F -.M^ 



VU{zifMV, 



dG 



VU{z2fMV, I ^G 



VC/(zi)^MV, ( ^ 

JK 



y VU{z,)^MW, { |G 



{xi,zi) 
{xl,zi) 

{X1,Z2) 
{xl,Z2) 



\ 



For a fixed n, we define the data 



/ ^^ \ 



.(1) 

Am 



By a least-squares method, we recover an estimation of the polarization tensor M{kn, D): 

M-':=argmin ||F(A/)-6„||. 

MeM 
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Again, once A'f(fc„,Z?) is estimated, a minimization approach yields correct parameter and size 
values. Since for any n, the eigenvectors of the matrix M{kn,D) are the ellipse axes, denoting r^ ^ 
and T^2 the estimated complex eigenvalues of AI(kn, D), one minimizes the following quadratic 



misfit functional 



E 

l<n<Af 



ab{kn- l){a + b) 



„est 



akn + b 



ab{kn - l){a + b) 



bkn 



^est 

'n,2 



over a, b, k, and e, in order to reconstruct the semi-axis lengths a and b and the material parameters 
k and e of the ellipse-shaped target D. 

If A^ is large enough, then semi-analytical formulas to estimate the semi-axis lengths a, b and 
the material parameters k, e hold. Since 



tIP} «7ra(a-|-6), 

one can estimate a and b as follows: 

^est 



^cst 



a 



'N,l 



^est 

'N,2 



^cst ^ 



7r6(a -I- 6), 



^est 

'N,2 



n T 



est 

AT.l 



^est 

'N,2 



\f^ 



^est 



(4.7) 



Table 4.2 gives estimations of a and b. The target is centered at zi — 1.5(cos(7r/3),sin(7r/3)) and 



the fish moves in the horizontal axis so that Z2 = (1.5cos(7r/3) — 1, 1.5sin(7r/3)). The material pa- 
rameters of the target are fc = 2 and e = 1. The data is collected for 10 frequencies equidistributed 
from 1 to 10. The reconstructed results are accurate. 
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Table 4.2: Estimations of the semi-axis lengths of ellipse-shaped targets using (4.7). 



Moreover, once the geometric parameters a and b are estimated, it is straightforward to recover 
k and e. Introduce 
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(1) 



^est 

'n,i 



f^n 1 



T:ab{a + b) a + knb 



From 



1 



isnujQ 



- a^i-n' 



,(1)' 



1 - b^J}r, 

one can estimate k and e as the real and imaginary parts of fc„. However, as shown in Figure [^9} 
one can see that the error on the real part is growing with the frequency. Therefore, in order 
to increase the robustness of the material parameter estimations, one estimate k using the lowest 
frequencies (for example the first three) and e using all the frequencies: 



1 ^ 



5R 
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a^'^Vn ^ ~ 



1 - 6^^*/il'' 



1 ^ 1 
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1 



-~J 



a°"*/ii'^ ~ 



UJqTI 



1 - 6°«V« ^ 



(4.8) 



Table 4.3 gives the material estimations using formula (4.8) for a disk and an ellipse. Once 



again the results are accurate. 
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Figure 4.9: Real and imaginary parts (respectively represented by squares and circles) for a disk- 
shaped target as functions of the frequency. Here, the target is with material parameters k — 2 
and e = 1, radius 0, 05, and placed at zi — 1.5(cos(7r/3), sin(7r/3)) and then at Z2 — (1.5cos(7r/3) — 
1, 1.5sin(7r/3)). The solid lines are the theoretical values. 
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Table 4.3: Estimations of the material parameters based on formula ( |4.8[ ). The disk has radius 
0, 05 and the ellipse has semi-axis lengths 0, 025 and 0, 1 and orientation angle 7r/3. Both targets 
are placed at zi = 1.5(cos(7r/3),sin(7r/3)) and then at Z2 = (1.5cos(7r/3) — 1, 1.5sin(7r/3)), and are 
illuminated with 10 frequencies equidistributed from 1 to 10. 



5 Conclusion 

In this paper, we have proposed a complex conductivity model problem for the quantitative anal- 
ysis of active electro-location in weakly electric fish. We have rigorously derived the boundary 
conditions to be used. We have proposed a non-iterative location search algorithm based on multi- 
frequency measurements. We have presented some numerical results which are promising. We 
have seen that increasing the number of frequencies (with not necessary different values) improves 
the stability. In fact, using multi- frequency measurements increases the signal-to-noise ratio. On 
the other hand, using different frequencies yields a faster robust location algorithm than repeating 
the data acquisition procedure with the same frequency. We have also proposed a procedure to 
reconstruct the electromagnetic parameters and the size of disk- and ellipse-shaped targets. This 
has been possible only because of multi-frequency measurements corresponding here to different 
frequency values. The use of multi-frequency measurements is fundamental in the characteriza- 
tion procedure. It has been known that polarization tensor for real conductivities cannot separate 
the size from material properties of the target [7j. For arbitrary-shaped targets, many important 
questions remain. In particular, it would be interesting to know how much parameter and size 
information one can extract from its polarization tensors for different complex conductivities. It is 
also worth mentioning that limiting our asymptotic expansions with respect to the target size to 
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the first-order term (the dipole approximation) does not give us the shape of the target. Hence, in 
a forthcoming work we will investigate how much information can be acquired in the near field by 
approaching the fish next to the target and developing the asymptotic expansions with high-order 
generalized polarization tensors [^. We will also investigate the stability of the proposed algo- 
rithm with respect to random fluctuations in the background permittivity and propose an original 
cross-correlation technique in order to correct for the effect of random heterogeneities on target 
location. 
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